Ring artifact and Poisson noise attenuation via volumetric multiscale nonlocal collaborative filtering of spatially correlated noise

Streaks and Poissonian noise in tomography data are attenuated using a new two-stage multiscale volumetric filtering framework where the degradations are modeled as correlated noise. The procedure is fully automatic, offers improved feature preservation, and can be further combined with regularized reconstructions to deliver state-of-the-art imaging quality.


Introduction
Computed tomography is commonly affected by streak noise in measured raw sinogram data (Jha et al., 2013;Artul, 2013;Boas & Fleischmann, 2012), which can be caused by miscalibration of detector linear response, beam fluctuations, beam hardening, or dusty or damaged scintillator screens (Haibel, 2008;Vidal et al., 2005;Anas et al., 2010). Streak noise in projections causes ring artifacts in reconstructed volumes, which present as centered circles or half-circles (Croton et al., 2019). As the sinogram data are obtained through a photoncounting detector, the statistics of the measured raw data can be further modeled through a Poisson distribution, which may result in high levels of Poissonian noise, commonly attenuated within the reconstruction process through iterative approaches (Mohan et al., 2014;Venkatakrishnan et al., 2013).
Although ring artifacts can be reduced by scanning protocols (Pelt & Parkinson, 2018), high-quality scintillator screens and detectors, it is difficult to completely avoid them and therefore achieve highest-quality reconstruction solely by experimental measures, requiring algorithmic processing of the acquisitions.
Popular methods to reduce ring artifacts include wavelet-FFT filters (Mü nch et al., 2009), combinations of polynomial smoothing filters and careful calibration of the detector response function (Vo et al., 2018;Croton et al., 2019), smoothing filters with segmentation in the tomogram domain (Massimi et al., 2018), ring removal in the tomogram domain upon polar coordinates transformation (Sijbers & Postnov, 2004;Li et al., 2021), and iterative algorithms (Paleo & Mirone, 2015) that combine regularized reconstruction with denoising. Recently, in Mä kinen et al. (2021), we proposed effective ring artifact attenuation through sinogram-domain collaborative filtering, presenting a multiscale architecture with a Block-matching and 3-D filtering (BM3D) image denoiser for correlated noise (Dabov et al., 2008;Mä kinen et al., 2020) at the core of the process. To the best of our knowledge, Mä kinen et al. (2021) offers state-of-the-art results in ring attenuation. In particular, it does not cause new artifacts around strong signal features, common to other popular ring removal algorithms. However, being based on a filter for 2-D data 1 , applied to individual sinograms, it may cause discontinuities across the third dimension.
In this work, we address both streak reduction and Poissonian noise removal from volumetric stacks of projections. The contribution of this work is threefold: (1) We propose a multiscale streak denoising framework for denoising of volumetric data. In particular, this framework can be seen as an extension of Mä kinen et al. (2021) for filtering of 3-D volumes.
(2) After streak noise removal, and before reconstruction, we embed a distinct multiscale denoising step to attenuate the Poissonian noise component of the projections. This allows to apply the reconstruction process using milder regularization and improve the tradeoff between noise reduction and artifact suppression.
(3) As a general-purpose algorithmic contribution, the filter used at the core of the multiscale denoising process is an improved version of the BM4D (Maggioni et al., 2012) volumetric denoising algorithm. The included enhancements, discussed in Appendix A, allow the long-range noise correlation which characterizes the streaks to be dealt with.
The proposed filtering procedure for both streaks and Poissonian noise is fully automatic and includes self-calibration of the filtering strength. We demonstrate the denoising performance of the proposed approach on real data from the table-top Prisma XRM microCT at Sigray, and from the synchrotron-based microCT at the Advanced Photon Source (APS) in Argonne, available through Tomobank (De Carlo et al., 2018).

Bright-field normalization
The following normalization of the raw projections and the streak model upon a logarithmic transformation follow that of Mä kinen et al. (2021).
The optical attenuation through the sample is determined experimentally via bright-field corrections through two separate acquisitions, the bright-field and the dark-field (Seibert et al., 1998). The bright-field is obtained by the imaging procedure with no sample, and the dark-field is obtained with no beam; both are 2-D arrays the size of effective pixels of the detector. The Beer-Lambert law further relates the X-ray transform through the sample to the optical attenuation by a logarithmic transformation (Swinehart, 1962).
Hence, the raw projections P raw are first normalized as where I D is the dark-field and I B is the bright-field 2 , and then log-transformed as 2.1. Noise model for normalized projections Apart from possible completely defective detectors 3 we treat the variation in detector response as normally distributed; as such, the streak noise will follow a normal distribution. Furthermore, we model the streak noise as locally stationary, meaning that the variance is presumed constant within the support of the denoising filter. Note that this does not mean that the noise is i.i.d. or white, as it is instead characterized by very long range correlation presenting as streaks.
As the data are obtained through a photon-counting detector, the statistics of the measured raw data can be further modeled through a Poisson distribution with nonstationary parameters after the bright-fielding.
Given the premises of normally distributed streak noise and Poissonian noise, noise in projections normalized by (1) can be formally written as where A are the underlying noise-free projections, P is the normally distributed streak noise component, and is Poissonian noise which we model as white and zero-mean; all components of (3) are considered as 3-D arrays and multiplications are elementwise.
The natural logarithm of (2) acts as a variance-stabilizing transformation (VST) for the multiplicative noise component ð1 þ P Þ. Hence, we have where the approximation comes from lnð1 þ P Þ % P .

Correlated noise
The denoising is conducted in two steps. First, we aim to estimate the 'streak-free' projections which are corrupted by white Poissonian noise. Then, as a separate denoising step, we consider the attenuation of the remaining noise originating from . Throughout this work, we will represent the volume to be filtered according to the correlated noise model presented in the following subsection. This model will assume different meaning at different parts of the algorithm. First, applying locally to the streaks as a type of long-range correlated noise; second, to the noise arising from the Poissonian component .

Correlated noise model
We consider the noisy input z : X ! R to be a combination of underlying data y and additive stationary spatially correlated noise to be filtered, where x 2 X & Z 3 is the coordinate in the finite three-dimensional volumetric domain X and being zero-mean i.i.d. Gaussian noise with unit variance, and Ã denoting 3-D convolution with the kernel g. The kernel g defines the spatial correlation of the noise as well as the noise strength, with kgk 2 = stdfg. An equivalent way of representing correlated noise is by its power spectral density (PSD) É, with F being the 3-D Fourier transform, and |X | denoting the cardinality (i.e. number of elements) of X. Equivalently, a kernel g satisfying (7)-(8) can be defined from É as

Estimation of noise standard deviation
When applying the above model to noisy data, it is essential to have knowledge of either the kernel g or, equivalently, the PSD É, as they fully characterize the noise. Assuming g in (7) is known modulo a scaling factor & from a known kernel g s , the noise estimation can be simplified to estimating &. In particular, in order to model the streak and Poissonian noise components arising from the particular composition of noise given in (3), the kernels g s should induce either very long range correlation or near white noise across each dimension d. The estimation procedure can be built as a direct extension of the one adopted by Mä kinen et al. (2021) to 3-D. To reduce the signal-to-noise ratio (SNR) to acquire a better noise estimate, we convolve z with a 3-D anisotropic kernel g d that provides either low-pass or high-pass filtering along different dimensions; g d is designed based on the noise statistics so that it preserves the noise component of interest while attenuating signal contrast. Specific instances are given in Section 4.1.1 and Section 5.3. One can then compute an estimate of the standard deviation of Ã g d via its median absolute deviation (Hampel, 1974), where smed denotes the sample median and the factor 1.4826 calibrates the estimate with respect to a normal distribution of the noise. As stdf Ã g d g = k&g s Ã g d k 2 , an estimate& & of & can be obtained througĥ

Multiscale streak filtering
In the following, we treat the first dimension of the stack of projections as the angular dimension, and the second and third as the horizontal and vertical displacement dimensions. Because the streaks are inherently low-frequency with respect to the angle, they are filtered entirely at a coarse angular scale; for this task, we extend the multiscale procedure of Mä kinen et al. (2021). The main changes in the proposed procedure arise from replacing the one-dimensional binning operations along the displacement dimension with corresponding 2-D binning operators B 2D and B À1 2D executed across both displacement dimensions. Furthermore, instead of using a direct 3-D extension of the 2-D streak PSD, we adjust the streak model to account for possible long correlation also along the displacement dimensions.
In detail, the multiscale streak attenuation procedure proceeds as follows. We begin by an angular binning B . The result of the angular binning Z 0 = B ðZÞ is binned K times through B 2D to obtain Z K = B K 2D ðZ 0 Þ. The size of each binned volume is a quarter of the input size. Then, we process each scale in a coarse-to-fine fashion, where progressively for each k = K; K À 1; . . . ; 2; 1; 0, we obtain an estimateŶ We start by taking as noisy input Z Ã K of BM4D the smallest binned volume Z K ; in this way, we obtain from Z Ã K = Z K the coarsest estimateŶ Y K , which is taken as initialization for the following recursive steps executed for each scale k = KÀ1; . . . ; 0: (1) Replace the coarser-scale components of Z k by those of the estimateŶ Y kþ1 : (2) Denoise Z Ã k with BM4D to produce the estimateŶ Y k . Finally, we replace the coarse angular components of the fullsize stack Z with those from the finest scale estimateŶ Y 0 ,

Multiscale noise model
For BM4D denoising, we regard Z Ã k of each scale k as z of the model (6), as where research papers and k = B k 2D ð 0 Þ = B k 2D ½B ð P Þ. This definition for Ã k , k < K, follows from considering the coarser-scale estimateŶ Y kþ1 as perfectly denoised. Similar to (8), Ã k is treated as correlated noise with PSD, where g Ã k is a correlation kernel and |X k | is the pixel size of Z k . As per (9), the kernel g Ã k can be defined as 4.1.1. Adaptive parametric model of É Ã k . We note that, in addition to the approximately white streak noise, the sinograms may contain streaks with very long range correlation across the displacement dimensions. As this correlation is aligned along the detector axis, it is not clearly observable in individual sinograms, but may create significant noise structure in the full volume. Hence, we approximate the streak noise 0 through three angularly constant streak noise components distinct in the displacement, where 0, w is streak noise white across both displacement dimensions, 0, u is streak noise constant across horizontal displacement, and 0, v is streak noise constant across vertical displacement. Let us denote by 0;p 2 f 0;w ; 0;u ; 0;v g a noise component of 0 . For each 0, p , we can define a respective scaled correlation kernel where & 0;p = stdf 0;p g, and kg 0;p k 2 = 1. Example realizations, kernels, and PSDs for each of these components as well as k are shown in Fig. 1 (top). We note that each k;p = B k 2D ð 0;p Þ is characterized by a kernel & k;p g k;p = 2 k & 0;p g 0;p . In particular, this property arises from the noise structure of the corresponding components: as B 2D operates through summation and elimination of adjacent pixels, the operation preserves both noise whiteness and constant noise. The ratio 2 k follows from the summation along two dimensions, meaning that the variance of the coarser scale is four times that of the finer scale.
Disregarding the specific support size of their actual finite realizations, we can identify the stationary random fields as where G, p is noise characterized by g k, p , and hence varf G;p g = 1. We can then express the residuals of any of the components 0, p as Then, where the noise correlation kernel corresponding to a where g B 2D is a 2-D kernel across the displacement dimensions characterizing the residual from 2-D binning of white noise. Specifically, where G, 2D is a two-dimensional white random field. The field size |X G | is included only for the normalization of the Fourier transform, and the formula holds for an arbitrary size. Top: example noise 0, p , the corresponding kernels g 0, p and the root PSD jF ½g 0;p j for each noise component in (17) with & 0;w = 6, & 0;u = 5, and & 0;v = 8, as well as example noise, kernel, and root PSD corresponding to the compound noise k . Bottom: example noise, the corresponding kernels, and root PSD of the corresponding binning residuals Ã k;p . For all visualizations, the angular dimension of the data is the vertical dimension in the figure. The DC corner of the Fourier spectra is marked by a circle. Note that all root PSDs are nonzero only on the angular DC plane, and the kernels and the noise consist of repeated planes across the angle.
Then, the PSD of Ã k , K ! k ! 0, can be written as As any Ã k is constant in angle, É Ã k is non-zero only across the DC plane with respect to the angular dimension. Example realizations, kernels, and PSDs for the residual components are shown in Fig. 1 (bottom).
Although (23) allows for modeling of very long range correlation, the streak noise is likely to contain minor correlation along the displacement not accounted for by this model. To adapt to such deviations, we allow the scaling parameters & k;p ! 0 for each noise component to vary with each scale k by estimating them individually at each scale, effectively accounting for mild local correlation in the noise.
Estimation of & k, w , & k, u , and & k, v . Based on (22) and (23), the PSD is completely determined by the values assumed by the three parameters & k, w , & k, u , and & k, v and the known kernels g k, p and g B 2D Ã g k;p . To adaptively obtain the parameters, we begin by obtaining three noise variance estimateŝ 2 k;w , 2 k;u , and 2 k;v . For each estimate, we define a corresponding filtering kernel g ðpÞ d such that 2 k;w estimates the variance of high-frequency streaks, 2 k;u estimates the variance of horizontally low-frequency streaks, and 2 k;v of vertically low-frequency streaks. For this purpose, we define d as a 1-D Gaussian function along dimension d, and d as a 1-D highpass kernel with Daubechies wavelet 'db3' of length 6 along d. Hence, convolution with d realizes low-pass filtering, and d realizes a high-pass filter. Then, g ðpÞ d is realized as a tensor product of three one-dimensional kernels across the dimensions d chosen based on the noise statistics through that dimension, Specifically, with m 0 , m 1 , m 2 as the pixel sizes of the three dimensions of Z Ã k , 0 is a 1-D Gaussian function along the angular dimension with standard deviation m 0 /8, and 1 and 2 are 1-D Gaussian functions along the two displacement dimensions with standard deviations of m 1 /12 and m 2 /12, respectively. Through these kernels, we obtain estimates of the three coefficients k;p as described in (10) and (11) with g s as either g k, p (k ¼ K) or g B 2D Ã g k;p (k < K).
We note that these three components do not directly correspond to & k, w , & k, u , and & k, v , as the frequencies of the white streak component Ã k;w partly overlap with those of Ã k;u and Ã k;v , i.e. Ã k;w includes also some low-frequency streak components. In particular, we have To this end, we can formulate a simple non-negative least-squares optimization as Finally, we construct the PSD through (22) and (23) with Variance of the streak noise may differ across the sinogram due to changes in photon flux or noise in the bright-field. Thus, it may not be possible to denoise Z Ã k assuming an equal É Ã k for all spatial positions without either oversmoothing or leaving noise artifacts in some areas. To adapt to nonstationarity, we further relax the streak model allowing the PSD to vary within each scale k. In particular, before noise estimation and denoising, we split Z Ã k into overlapping, volumetric segments. We apply BM4D separately on each segment of Z Ã k , using a PSD scaled by parameters estimated from the same segment, i.e. we consider each segment as a separate noisy volume z with a corresponding É. After denoising, the segment estimates produced by BM4D are recombined with a windowing function to form the full estimateŶ Y k .

Attenuation of extreme streaks
We note that the projections often include several streaks caused by defects in the scintillator. These streaks can be far stronger than that reasonably produced by the distribution of P and therefore require a specific pre-processing. To this end, after the bright-fielding and log-transform and before the multiscale denoising procedure, we apply the simple extreme streak attenuation procedure as described by Mä kinen et al.
(2021), which applies median filtering on extreme streak values detected through local polynomial fit of angular medians.

Poisson denoising
A filter for additive noise is not immediately applicable to the approximately white noise ofŶ Y originating from the Poissonian component . Firstly, the bright-fielding (1) introduces substantial spatial variability in the Poisson model. As a result, for a given optical attenuation, noise in bright-fielded projections can be stronger or weaker in different parts of the detector, for example around edges in cone-beam acquisition. Secondly, while the logarithm (2) effectively makes the streak noise additive, it also changes the typical affine-variance model of the Poissonian noise to a nonlinear one where the variance is not constant, but asymptotically inversely proportional to the mean. In order to model the noise inŶ Y through (6), we take care of these two issues as follows.

Reducing nonstationarity induced by bright-fielding
The Poissonian noise component originates from a counting process which takes place before bright-fielding (1), and specifically before the division by I B ÀI D , which intro-research papers duces a spatially variant scaling of the variances. To undo this scaling, we consider where I L = I B À I D . Then, S can be treated as the log-scale version of a homogeneous Poissonian process; S is thus subject to signal-dependent noise where the variance of the noise can be expressed as a smooth nonnegative function of the underlying signal, where the same F applies to each pixel. In particular, it can be shown that asymptotically for large flux F EfSg ð Þ / EfSg ! 1 1 EfSg :

Stabilization of variance
To turn a model like (25) into (6) we again resort to the use of a VST. As large-flux asymptotics are irrelevant for denoising problems characterized by low signal-to-noise ratio, and to pragmatically accommodate for model uncertainties, we model F as a polynomial with arbitrary data-driven coefficients. The method (Foi, 2009;Azzari & Foi, 2014) simultaneously identifies the coefficients for an arbitrary signaldependent noise model where the variance is a positive power of an unknown polynomial, and returns the associated variance stabilizing transformation f as well as the corresponding exact-unbiased inverse VST f À1 EUI . An example of an estimated standard deviation function ðF FÞ 1=2 and the corresponding VST f are illustrated in Fig. 2, where the effectiveness of the stabilization can be deduced by the estimates of d stdf f ðSÞg stdf f ðSÞg being scattered around 1.

Multiscale denoising of the stabilized Poisson noise
To avoid introducing structured artifacts that are present in the bright-field and dark-field images, we further consider a smoothed version I smooth L of the field component, where g I is a 2-D Gaussian kernel, and medfilt denotes a 2-D median filter. The median filter is adopted in order to remove extreme outliers (e.g. from broken pixels), and the convolution with the Gaussian ensures a smooth result. Then, I smooth L can be used for approximate correction for the brightfield induced nonstationarity witĥ Y Y þ lnðI smooth L Þ. The stabilized noisy stack can then be written as where corresponds to the stabilized noise andỸ Y to the signal upon stabilization.
We consider white, and assume the streak denoising procedure to remove all streak noise frequencies, including those of . Hence, we treat as missing the streak frequencies, i.e. with a PSD, where c is a constant such that varf g = cjXj À1 ðm 0 À 1Þ m À1 0 . For multiscale denoising of the Poisson component, we define three-dimensional binning and debinning functions as B 3D = B 2D B and B À1 3D = B À1 B À1 2D , and obtain K Poi scales of binned noisy volumes asZ Z k = B 3D ðZ ZÞ, k 2 f0; . . . ; K Poi g. Then, unlike the progressive denoising of the streaks, we begin by BM4D denoising ofZ Z k of each scale k; at each scale, we model the noise through a PSD of the form (28). This way, we obtain an initial estimateŶ Ỹ Y Y k of the corresponding noise-free volumẽ Y Y k at each scale. Then, starting from k ¼ K À 1, we combine only the denoised volumes of each scale by recursively replacing the low-scale components ofŶ Ỹ Y Y k , k ¼ fK À 1; . . . ; 0g, by those of the lower scale, where g G denotes a 3-D Gaussian kernel. Although the low frequencies are obviously denoised more effectively in the coarser scale, the higher frequencies of the coarser scale are commonly estimated worse than the respective estimate of the finer scale (Facciolo et al., 2017). As such, g G realizes a lowpass filtering which selects only low frequencies of the coarser estimate to be used in the full estimate.
To account for possible remaining nonstationarity and slight correlation of the noise, we perform the denoising in segments similar to as described in Section 4.1.2 for streak noise, and estimate a separate scaling parameter & 2 Z Z in construction of the PSD at each scale. In particular, we estimate & 2 Z Z as described in Section 3.2 with g d = 0 1 2 and g s = As (30) negates the field correction, we note that had we used the non-smooth field I L in (27) [and respectively in (30)], any noise or spurious structures present in I L could be introduced into d lnðAÞ lnðAÞ, as they might have been denoised by BM4D and hence not preserved inŶ Ỹ Y Y Ã 0 . Upon variance stabilization, Poissonian data become asymptotically normal (Curtiss, 1943). Due to the additional Gaussianization induced by the binning and by the linear transformations operated by the filter, the assumption of normality in (7) can be adopted for denoising of the Poissonian component in this work even for low-count data.
The full denoising process is shown in Fig. 3.

Experiments
We test our pipeline on synthetic data as well as two real acquisitions displaying ring artifacts and Poisson noise.  (Cocosco et al., 1997) obtained through a padding and Radon transform upon a sign change and an exponential transformation. Specifically, we regard this stack as the underlying projections A and generate noise according to (3) with g as a constant of size m 0 Â1Â1 (equal to g 0, w of Fig. 1). To obtain streak noise of different strengths, the streak noise component ð1þ P Þ is generated with stdf P g = 0.005, 0.01, 0.02, 0.05. Next, to generate noisy measurements with different SNR levels for the Poisson component, we separately scale A to the ranges [2560, 5120] (higher SNR), [1280,2560], and [640, 1280] (lower SNR) and generate a Poisson variate with mean and variance Að1þ P Þ, thus defining the Poissonian noise as the difference between this Poisson variate and Að1þ P Þ. Furthermore, we include experiments with = 0 (infinite SNR), thus resulting in a total of 16 combinations of Poisson and streak noise strengths. We do not simulate extreme streaks or the bright-fielding (I B = 1 and I D = 0). For the streak removal, we consider ln½A þ =ð1þ P Þ as the streak-free yet noisy stack Y.
The results of the phantom experiments 4 for streak attenuation are collected in Table 1, and, for full denoising, evaluating the reconstructed volumes, in Table 2 using iterative regularized TV reconstruction with optimized regularization parameter strength r. The experiments for both streak and Poisson denoising are illustrated in Figs. 4 and 5. All reconstructions are performed upon a sign change.
The Fly dataset consists of 180 projections with 50 s exposure (detector pixel size 27 mm, demagnified to 15.7 mm by cone-beam geometry) collected using a Sigray Prisma X-ray micro-tomography instrument at 34 kV; the detector size is 512Â512 pixels. The denoising results for two different sinograms, as well as a corresponding tomogram after streak The full denoising process, requiring as inputs the noisy projections P raw and the bright-and dark-fields I B , I D , (1) and producing as the output an estimate d lnðAÞ lnðAÞ (30) of the underlying stack of projections ln(A) (3). As an intermediate output, an estimateŶ Y (12) of the streak-free yet noisy stack of projections Y (5) is also produced. attenuation, are shown in Fig. 6. A comparison of denoising on a vertical slice of the stack of tomograms is shown in Fig. 7, and a comparison for fully denoised reconstructions is shown in Fig. 8.
We also test the algorithm on a soft tissue sample 00072 displaying severe ring artifacts freely available in TomoBank (De Carlo et al., 2018). The data contain 1500 projections with 1.43 mm pixels, obtained at the Advanced Photon Source, 2-BM beamline; other experimental parameters are X-ray energy of 20 keV, 10 mm LuAG scintillator, and sampleto-detector distance of 15 mm. The detector size is 2160Â2560 pixels. Included are ten samples for bright-and dark-fields, which are averaged to obtain a single bright-field and dark-field. A sinogram and a corresponding tomogram from the denoising results for streak removal are shown in Fig. 9, and slices of the stack of tomograms are compared in Fig. 10. Reconstructions of fully denoised projections are further compared in Fig. 11.  Table 2 Average SNR for the reconstructed volumes of the BrainWeb phantom for the set of experiments shown in Table 1.   Table 1 Average signal-to-noise ratio for attenuation of streaks in the BrainWeb phantom subject to mixed streak and Poissonian noise as in (3), with different combinations of stdf P g and peak values of A, with 'peak' = 1 being the limiting case for which = 0.
Left-to-right: noisy stack of projections Z (4), and estimatesŶ Y of the stacks of projections Y (5) denoised by the proposed procedure (12) Vo et al. (2018). As the table compares only streak removal, the SNR values are calculated with respect to the streak-free yet noisy projections Y = ln½A þ =ð1þ P Þ as SNRðŶ YÞ = 10 log 10 fsvar X fYg= smean X ½ðŶ Y ÀYÞ 2 g, where svar and smean denote sample variance and sample mean, respectively. Each value of the table is the average SNR over 10 different noise realizations.
Munch et al.
Vo et al. The proposed method achieves superior SNR values in streak removal in all simulated noise experiments. Although the difference to the 2-D implementation of Mä kinen et al.  Vo et al. (2018), the proposed method avoids creation of shadow artifacts around strong signal features. Furthermore, performing the Poisson denoising through the proposed framework allows application of standard filtered back-projection reconstruction to data originally corrupted by Poisson noise, but can also improve quality of iterative TV reconstruction.

Parameters
For streak attenuation, we calculate K following the formula of horizontal binning from Mä kinen et al. (2021), using as the base the size of the smallest displacement dimension. As a result, we use K = 5 for 00072, K = 3 for Fly, and K = 2 for the phantom. These values were found to offer a reasonable compromise between denoising wide streaks versus preserving low-frequency signal components. Other processing parameters are adjusted for the smaller block size and processing neighborhood of BM4D. For angular binning, we use m = dm=dm=32ee % 32 pixels, where m is the original angular size and m the output size; the resulting size is half of that used by Mä kinen et al. (2021). For segmentation of the streak denoising, we use a window of size dm =2eÂ19Â19 pixels. For the Poisson denoising, we use K = 1 and mÂ19Â19 segments. For variance stabilization, we use the implementation ClipPoisGaus (Azzari & Foi, 2015) of Foi (2009) and Azzari & Foi (2014), and use a quadratic polynomial for the variance model F.

Discussion and conclusions
We have presented a framework for three-dimensional attenuation of streak noise extending the 2-D framework of Mä kinen et al. (2021), as well as a BM4D denoiser utilizing the algorithmic improvements of Mä kinen et al. (2020). Furthermore, we have included a denoising step for Poisson noise in the sinogram domain through BM4D and the adaptive variance stabilization of Foi (2009) and Azzari & Foi (2014).
We test the algorithm on both synthetic and real data, demonstrating superior SNR compared with other popular streak removal algorithms, and showing improvements in streak attenuation over Mä kinen et al. (2020). Furthermore, we compare the results with those which use the conventional BM4D for correlated noise, demonstrating that the included improvements in BM4D for correlated noise are essential for successful streak attenuation. The included Poisson denoising allows for full sinogram-domain denoising within the framework. By operating fully in the 3-D stack of projections, the 3-D structure of the data can be leveraged for more accurate noise removal. The proposed procedure is fully automatic and does not require extra input parameters.
To compare different methods under their own ideal conditions, we have specifically selected the TV regularization parameter values that provide the best reconstruction quality. However, in real-world applications, these values cannot be identified precisely, and too small or too large parameter values may lead to residual noise or excess smoothing of the reconstructions. Inclusion of the proposed Poisson denoising step allows for weaker regularization, but notably also reduces the effects of relative shifts in the parameter values, meaning that the reconstruction can be safely deployed even when the regularization cannot be precisely tuned.
To consider the computational cost, we note that both denoising steps of Fly (181Â512Â512 pixels) run single-  threaded on an AMD Ryzen 7 1700 processor each take about one hour. The computational cost is mostly due to the BM4D denoising in CPU. Although the adopted implementation is single-threaded, the algorithm is embarrassingly parallel, and thus a highly parallel GPU-based implementation is expected to reduce the total run time to the scale of seconds (Davy & Ehret, 2020).
The Poissonian noise attenuation can also be performed without the preceding ring reduction step on data which do not display ring artifacts. In such case, ÉZ Z should be replaced by a flat PSD, as the Poissonian noise is approximately white Comparison of reconstructions of the 3-D BrainWeb phantom corrupted with streak and Poisson noise as in (3), corresponding to the sinograms shown in Fig. 4. Top: ground truth volume, and reconstructions of ln(A) (3) and Z (4) obtained through filtered back-projection. Bottom: comparison of TV reconstruction of estimates with various regularization strengths r, where the percentage implies a multiplier to the regularization optimized to maximize SNR, i.e. 100% means 'SNR-optimal' regularization. Top-to-bottom: proposed full estimate ln(A), proposed streak-free estimateŶ Y, and streak-free estimate ofŶ Y (Vo et al., 2018), each with 100%, 50%, and 150% relative regularization strengths. Proposed estimates are computed embedding the improved BM4D. Notably, the full filtering offers improved reconstruction quality, and is also less sensitive to variations in the regularization parameters.

Figure 6
Denoising of the stack of projections of Fly, showing two sinograms of the noisy stack of projections Z (4) and the corresponding estimates d lnðAÞ lnðAÞ of the underlying stack of projections ln(A) (3) obtained with the proposed framework (top), and the tomograms of the second sinogram (bottom), obtained with filtered back-projection using cone-beam geometry (Feldkamp et al., 1984), for both the noisy data Z and the proposed estimateŶ Y of the streak-free projections Y (5). The tomogram for d lnðAÞ lnðAÞ is shown in Fig. 8. The first sinogram shows significant model nonstationarity in both streaks and the Poissonian component due to the bright-fielding.    prior to streak attenuation, whereas (28) considers the streak noise frequencies removed. Running the full denoising procedure in the absence of either streak or Poisson noise will lead to very small estimates for the corresponding noise components, meaning that no significant denoising will be performed for that noise.
We note that although we have focused on the full denoising of the projections, typical reconstruction pipelines, such as the iterative TV, provide further noise attentuation. For best results in combining the proposed denoising procedure with  Denoising of the stack of projections of 00072. Top-to-bottom: display of a single sinogram of the noisy stack of projections Z (4), the corresponding estimate of the underlying projections d lnðAÞ lnðAÞ from the proposed procedure, and the corresponding tomograms of Z and the estimatesŶ Y of streak-free stacks Y (5), respectively, obtained with filtered back-projection. The tomogram for d lnðAÞ lnðAÞ is shown in Fig. 11. Although the data present challenges through inconsistent noise intensities across the angular dimension, most streak noise and Poissonian noise is attenuated without notable loss of signal.  such pipelines, it may be necessary to adjust the filter strength for the denoising of Poissonian noise, e.g. for reduced attenuation of high-frequency noise, as it is further attenuated within the reconstruction. Likewise, integration of the proposed procedure within an iterative alternating reconstruction is left for future study.
APPENDIX A Collaborative filtering and the BM4D denoising algorithm A1. Collaborative filtering The rationale of transform-domain filtering is to work with a representation of the signal where most of the signal is compacted to only a few coefficients, whereas the remaining coefficients mostly comprise noise. Hence, by attenuating the coefficients with a non-linear shrinkage operator, it is possible to attenuate noise while keeping most of the signal intact. Nonlocal collaborative filters utilize this property in the context of collective transform coefficients of groups of similar patches extracted from the input. One of the most popular collaborative filters is the Block-Matching and 3-D filtering (BM3D) (Dabov et al., 2007) denoising algorithm, which performs denoising on groups of blocks extracted from a 2-D image. In the BM4D volumetric denoiser (Maggioni et al., 2012), the patches are 3-D volumes extracted from the volumetric data.
All operations of collaborative filters are made with regard to a reference patch moving through the volume. For each position of the reference patch, the following steps are executed: (1) Collect similar patches into a group through patchmatching.
(2) Obtain a group transform spectrum by collectively transforming the group of patches.
(4) Transform the shrunk spectra back to patch estimates and aggregate them to the original locations from which they were collected.
For details about the algorithm in arbitrary dimensionality, we refer the reader to Mä kinen et al. (2020). In the following section, we describe special considerations for the implementation of the algorithmic improvements in Mä kinen et al.  Comparison of fully denoised tomograms of 00072, corresponding to the sinogram in Fig. 9. Top-to-bottom: tomogram reconstructed from a noisy sinogram of the stack of projections Z (4), of the estimate of stack of underlying projections d lnðAÞ lnðAÞ with the proposed procedure with FBP reconstruction, and tomograms of the estimates for streak-free stacksŶ Y of Mü nch et al. (2009) with TV reconstruction, andŶ Y of Vo et al. (2018) with TV reconstruction. TV regularization was tuned visually, balancing residual noise and smoothing of signal. Compared with the reference methods, the proposed procedure manages to remove most streaks without significant loss of detail, as well as most Poissonian noise without excess smoothing of the signal. simple procedure which obtains a 2-D projection of the 3-D PSD by preserving the two largest principal components, aiming to preserve the characterizing features of the PSD shape. This projection is then used to compute features as described for a 2-D PSD by Mä kinen et al. (2020) for the estimation of suitable and 2 .
A2.2. Fast implementation. We consider the fast implementation as suggested by Mä kinen et al. (2020). In particular, we perform all operations on a downscaled PSD of size N f ÂN f ÂN f and compute exactly only the K f first volumes of the 4-D spectrum and approximate the rest using the conventional variances. Furthermore, Fourier symmetries and sparsity of the transformed arrays can be exploited to reduce computational cost similar to the 2-D case.
A2.3. Refiltering. As noted by Mä kinen et al. (2020), even with exact modeling of the collaborative transform-domain noise spectrum, the accuracy of collaborative filtering is limited by the systemic factors arising from the used transforms, both in size and possible symmetries of the transform spectrum which may limit the modeling of the global PSD. As a result, the denoising may attenuate excess signal, leading to oversmoothing in some frequencies; Mä kinen et al. (2020) proposes the mitigation of these systemic issues through an extra filtering step performed on the denoising residual. The three-dimensional spectra are not exempt from these limitations, and as such we adopt the global Fourier thresholding and refiltering procedure through a 3-D FFT.